!> final output prints
!!
!! this module has all data and procedures for defining and
!!  manipulating the problem output, namely the fluxes
!!  and eigenvalues
module class_output

    use class_input, only: indat, uout
    use class_dat, only: rhsdat,coefdat,flux,keff

    implicit none

    contains

    !> drives output processes
    subroutine output(inp,phi)
#include "finclude/petsc.h90"
        type(indat), intent(in)                  :: inp
        type(flux), intent(in)                   :: phi
    	! if we want to print fluxes
        if (inp%printout.eq.1) call printflux(inp,phi)

    	! if we want to make a silo file for printing
        if (inp%plotout.eq.1)call plotflux(inp,phi)

    end subroutine output

    !> print the fluxes
    subroutine printflux(inp,phi)
#include "finclude/petsc.h90"
        type(indat), intent(in)                  :: inp
        type(flux), intent(in)                   :: phi
        PetscInt :: ni, i, ierr
        PetscInt, pointer :: ix(:)
        PetscScalar, pointer :: y(:)
        write (uout,*) " preparing to print fluxes "
!PetscErrorCode  VecGetValues(Vec x,PetscInt ni,const PetscInt ix[],PetscScalar y[])
!Input Parameters
!	x 	- vector to get values from
!	ni 	- number of elements to get
!	ix 	- indices where to get them from (in global 1d numbering)
!Output Parameter
!    y  - array of values
        ni = inp%nxfm * inp%nyfm
        allocate( ix(ni), y(ni) )
        do i = 1, ni
            ix(i) = i-1
        end do
        do i = 1, inp%numg
            ! get the group flux
            call VecGetValues( phi%val(i), ni, ix, y, ierr )
            ! print the group flux
            write (uout,*) "  group: ",i
            write (uout,'(5e12.5)') y
        end do
        deallocate(ix,y)
    end subroutine printflux


    !> create silo file for plotting fluxes
    subroutine plotflux(inp,phi)
        include "silo.inc"
        type(indat), intent(in)                  :: inp
        type(flux), intent(in)                   :: phi
        PetscInt                                 :: ierr, err
        integer                                  :: dbfile, i
        PetscScalar                              :: x(inp%nxfm+1), y(inp%nyfm+1)
        PetscScalar, pointer                     :: dx(:), dy(:), f(:)
        character(len=20)                        :: tmp='____________________',tmp2='____________________'

#include "finclude/petsc.h90"

         ierr = dbcreate("fluxes.silo", 11, DB_CLOBBER, DB_LOCAL, &
                    "group flux data", 22, DB_PDB, dbfile)
         if( dbfile.eq.-1 ) then
             write (*,*) "Could not create Silo file!"
         end if

         ! pull dx and dy scalar values from the input vectors inp%dx and inp%dy
         call VecGetArrayF90(inp%dx,dx,ierr)
         call VecGetArrayF90(inp%dy,dy,ierr)

         ! from the scalar arrays dx and dy and the initial val's of xcm and ycm
         ! build the x and y meshes
         x(1) = inp%xcm(1)
         y(1) = inp%ycm(1)
         do i = 2, inp%nxfm+1
             x(i) = x(i-1)+dx(i-1)
         end do
         do i = 2, inp%nyfm+1
             y(i) = y(i-1)+dy(i-1)
         end do

         ! restore the dx and dy vectors
         call VecRestoreArrayF90(inp%dx,dx,ierr)
         call VecRestoreArrayF90(inp%dy,dy,ierr)

         ierr = dbputqm ( dbfile, "quadmesh", 8, "xc", 2, "yc", 2,               &
                          "zc",2, x, y, DB_F77NULL, (/inp%nxfm+1,inp%nyfm+1,1/), &
                          2, DB_DOUBLE, DB_COLLINEAR, DB_F77NULL, ierr           )

         do i = 1, inp%numg
             ! make a group-specific string
             write (tmp,*) i
             tmp2='group' // adjustl(tmp)
             write(*,*) tmp2

             ! get group flux
             call VecGetArrayF90(phi%val(i),f,ierr)
             ! write it to the database
             ierr = dbputqv1( dbfile, tmp2, len_trim(tmp2), "quadmesh", 8, f,                 &
                              (/inp%nxfm,inp%nyfm,1/), 2, DB_F77NULL, 0,         &
                              DB_DOUBLE, DB_ZONECENT, DB_F77NULL, err            )
             ! restore the flux
             call VecRestoreArrayF90(phi%val(i),f,ierr)
         end do
         ierr = dbclose(dbfile)
    end subroutine plotflux

end module class_output
